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Numerical  Approximations  In  Heat  Transfer  Problems  and 
Usage  of  IBM  7090  Computer  For  Solutions 

I.  Introduction 

The  purpose  of  this  paper  Is  to  furnish  Information  necessary  to 
obtain  numerical  approximations  to  the  solutions  of  certain  heat  transfer 
problems  by  making  use  of  a  thermal  model  and  related  IBM  7090  computer 
program  developed  by  BBE  and  BCG.  The  solutions  take  the  form  of  a  time 
history  of  temperature  distribution.  The  problems  concern  heated  structures 
or  components  of  quite  general  geometric  configuration  and  material  composi¬ 
tion. 

Material  properties  may  be  temperature  dependent.  Heating  may 
occur  by  convection,  by  radiation,  or  by  conduction.  In  the  case  of  aero¬ 
dynamic  heating,  provision  Is  made  for  a  real  gas  boundary  layer  with  auto¬ 
matic  determination  of  whether  the  flow  Is  laminar  or  turbulent.  The  tra¬ 
jectory  and  ambient  atmosphere  may  be  arbitrarily  specified. 

In  section  II,  we  consider  a  thermal  model  that  consists  of  a 
lumped  parameter  network  of  thermal  capacitances  and  conductances.  This 
model  may  be  thought  of  as  arising  from  the  replacement  of  the  differential 
equations  of  heat  conduction  and  their  boundary  conditions  by  a  suitable  set 
of  difference  equations  or  from  the  replacement  of  a  continuous  heat  conduc¬ 
tion  system  by  a  corresponding  lumped  parameter  system  on  the  basis  of  physical 
considerations.  Our  goal  here  Is  not  to  discuss  the  limitations  of  the  lumped 
parameter  model*,  but  rather  to  describe  the  equations  which  govern  it.  In 
the  appendices  we  give  first  the  details  of  the  computations  Involved  In  the 
calculation  of  aerodynamic  heating.  Next,  we  give  a  criterion  to  evaluate 
the  stability  of  the  calculations.  Then,  we  show  the  difference  equations. 

In  section  III,  we  consider  the  specification  of  this  model  via  a 
very  simple  but  complete  FORTRAN  control  program.  The  existence  of  a  certain 


*For  an  excellent  treatment  of  the  derivation  of  lumped  parameter  heat  transfer 
models  see  Dusenberre's  Numerical  Analysis  of  Heat  Flow,  McGraw-Hill  Book  Co., 
Inc. ,  1949. 
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set  of  FORTRAN-FAF  subroutines  In  the  BCC  Library  (Nos.  11.02.01-11.02.12)  \/ 

then  makes  direct  liaison  between  the  BCC  Operations  Project  and  the  user 
a  relatively  simple  matter  which  requires  virtually  no  knowledge  of  pro¬ 
gramming  on  the  part  of  the  user.  The  dimensional  units  to  be  used  In  each 
case  are  listed  In  Appendix  D. 

Given  a  specific  problem,  one  proceeds  by  setting  up  the  model 
along  the  lines  Indicated  In  section  II  and  specifying  the  model  along  the 
lines  Indicated  In  section  III.  In  section  IV  the  printout  format  1j  dis¬ 
cussed.  Finally,  a  t3rplcal  example  Is  carried  through  In  detail  In  section 
V.  The  last  appendix  gives  the  limitations  on  the  size  of  the  program  and  the 
estimates  which  are  to  be  used  In  computing  the  machine  time. 

II.  The  Thermal  Model  and  Its  Equations 

The  thermal  models  that  we  consider  consist  of  networks  of  thermal 

capacitances  and  conductances.  In  any  model,  to  each  end  point  of  a  thermal 

conductance  Is  assigned  a  number  from  1  to  1000  called  an  Index.  Each  such 

point  Is  called  an  Indexed  point  and  may  have  assigned  to  It  a  value  of 

thermal  capacitance  or  a  temperature  prescribed  as  a  function  of  time.  In 

any  case  the  temperatures  which  are  associated  with  the  Indexed  point  s,  for 

example,  are  designated  T  .  .  Similarly,  the  value  of  a  thermal  capacitance 

s 

associated  with  an  Indexed  point  r  Is  designated  C^,  and  the  value  of  a 

conductance  joining  the  Indexed  points  r  and  s  Is  designated  by 

general  the  values  of  each  thermal  capacitance,  C^,  will  be  dependent  on  the 

temperature,  T  ;  and  the  values  of  any  thermal  conductance  K  will  depend  on 
S  ITS 

T^,T^,  and  time.  By  allowing  this  generality  the  thermal  capacitances  can 
represent,  over  an  extended  temperature  range,  the  thermal  capacities  of 
pieces  of  solids;  and  the  thermal  conductances  can  represent  the  behavior  of 
the  thermal  conductivities  of  solids,  of  convective  heat  transfer,  and  of 
thermal  radiation. 

Within  the  network  at  any  time  there  will  be  a  thermal  flux  between 

each  pair  of  points  joined  by  a  thermal  conductance.  The  value  of  the  flux, 

q  ,  from  r  to  s,  for  example.  Is  given  by  an  apparently  linear  relation, 
rs 


•rs 


rs 


(T  -  T  ) 

I  s 


(1) 
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The  linearity  of  this  relation  is  truly  only  apparent,  since  may  be 

temperature  dependent.  To  determine  the  total  thermal  flux,q  ,into  the 

s 

point  s,  It  is  necessary  only  to  sum  all  the  q's  which  have  s  for  their 
second  index.  This  is  expressed  by  writing 


r 


Further,  in  a  short  period  of  time  At,  the  total  thermal  energy 

into  the  point  s  is  q  .  At*  Consequently,  if  s  is  a  point  with  which 

8 

a  thermal  capacitance  C  is  associated,  then  the  change  in  temperature  AT 

s  s 

in  Interval  At  is  computed  by  the  relation 

q  At 

■  -V  ■ 

s 

Here  C  is  presumed  to  depend  on  the  value  of  T  at  the  beginning  of  the 
s  s 

time  Interval,  the  precise  nature  of  this  dependence  will  be  discussed  later. 

If  the  temperatures  of  all  Indexed  points  are  known  at  a  time  t 

and  if  the  values  of  the  C  and  K  are  also  known,  then  the  relations  (1), 

s  s 

(2) ,  and  (3)  give  the  values  of  the  temperatures  of  all  indexed  points  which 
have  capacitances  associated  with  them  at  a  slightly  later  time  (t  +  At). 

At  the  remainder  of  the  Indexed  points  the  temperatures  at  (t  +  At)  are 
prescribed  or  calculated  directly  as  functions  of  time,  so  that  by  step  by 
step  computations  the  temperatures  can  be  obtained  at  all  Indexed  points  at 
later  times. 

The  general  description  of  the  thermal  network  must  be  completed 
by  prescribing  the  computation  of  the  values  of  the  conductances  and  capa¬ 
citances.  Of  these  elements  the  simplest  in  form  is  the  thermal  capacitance. 
For  any  of  these,  say  C^,  the  defining  equation  is 

“s  ■  '’s  •  V.  ■ 


where  V  is  a  prescribed  constant  having  the  dimensions  of  volume,  and  (PC  ) 

®  P  s 

is  given  either  by  a  polynomial  in  degree  five  or  less  or  by  Interpolating 


in  a  table,  and  has  the  dimension  of  thermal  capacity  per  unit  volume. 
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The  thermal  conductances  may  be  of  any  of  several  types.  Of 
these  the  simplest  has  the  form 


K 


rs 


rs 


(5) 


where  A  and  L  are  the  area  and  length  respectively,  and  Is  a  thermal 
conductivity  which  depends  on  tenqserature.  In  particular,  k  Is  given  by 

ITS 

a  polynomial  of  degree  five  or  less  or  by  Interpolation  In  a  table  and  Is 
evaluated  at  the  average  temperature,  \  (T^+T^). 

A  more  complicated  thermal  conductance  Is  a  composite  of  two  of 
the  simple  types  just  discussed  joined  by  a  "contact"  conductance.  Here  It 
Is  convenient  to  Introduce  two  additional  network  points  (suppose  r'  and  s') 
at  the  ends  of  the  contact  resistance  as  Is  shown  In  the  following  diagram: 


• - V\A/ — - ^AAy - - - VV\/ — • 

r  r'  s'  s 

It  should  be  noted  here  that  the  points  r'  and  s'  are  not  regarded 
as  Indexed  points  since  they  are  Internal  to  a  single  conductance.  Since 
there  are  no  lumped  thermal  capacitances  associated  with  the  points  Indexed 
by  r'  and  s',  the  thermal  flux  through  each  conductance  must  be  the  same  as 
the  flux  from  r  to  s.  Thus  the  relations  are 


q 

q 

q 

q 

q 


rs 

'^rr'  **r's' 

“  %'s 

rs 

- 

K  (T  -  T  )  , 
rs  r  s 

rr' 

- 

K  ,  (T  -  T  ,) 
rr'  r  r' 

> 

r's 

1  * 

r  8  r'  s 

.)  , 

CD 

CD 

- 

K  ,  (T  ,  -  T  ) 
s's  s'  s 

• 
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Together  these  lead  to  the  formulation: 

1 


rs 


K  , 
rr 


K  ,  ,  K  , 
r  8 '  s  s 


(6) 


T 


r' 


T 

s 


•rs 

K  , 
s  s 


(7) 


Here  the  contact  conductance  the  same  form  as  that  given  by 

(5)  except  that  A/L  is  replaced  by  A  alone.  As  suggested  by  the 

diagram  K  ,  is  calculated  from  the  mean  of  T  and  T  , ,K  ,  ,  from  that 
o  rr'  r  r'  r  s' 

of  T  ,  and  T  , ,  and  K  ,  from  that  of  T_ ■  and  T 
r'  s'  ’  s's  8  8* 

In  order  to  compute  the  value  of  a  composite  conductance 
the  values  of  its  component  conductances  must  be  'xnovm.  These  in  turn 
depend  on  the  internal  temperatures  T^,  and  T^,  as  well  as  on  T^  and  T^. 
However,  the  values  of  T  ,  and  T  ,  depend  on  T  and  T  and  the  values  of 
the  conductances.  This  poses  an  iiq>licit  non-linear  problem  for  the  com¬ 
putation  of  T  ,,  T  ,,  and  the  value  of  the  composite  conductance.  To 
IT  8 

solve  this  problem  with  adequate  accuracy  at  each  time  step  the  following 
iterative  procedure  is  used.  The  values  of  T  ,  and  T  ,  at  the  beginning 

IT  8 

of  the  last  time  step  and  the  just  calculated  values  of  T^  and  T^  at  the 
end  of  the  time  step  are  used  to  compute  an  approximating  value  for  the 
component  conductances  and  thus  for  the  composite  conductance.  From  this 
value  of  the  composite  conductance,  T^,  and  T^  a  value  of  ,  and  thus 
new  values  of  T^,  and  T^,  can  be  computed.  These  new  values  are  then  used 
to  recompute  the  values  of  the  conductances.  This  procedure  may  be  repeated 
as  many  times  as  necessary  in  order  to  obtain  good  values  of  the  composite 
conductance  and  the  internal  temperatures. 

A  third  type  of  thermal  conductance  that  we  consider  arises  from 

thermal  radiation.  In  this  case  q  is  given  by 

rs 


q 


rs 


CTAF.  F,  (T^  -  tS 
As  '  r  s 


(8) 
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where , 


o  “  Stefan-Boltztnann  constant 
A  *  surface  area  of  element  in  question 


F*  (T  .T  ) 


In  this  case  we 


geometric  exchange  factor 

net  emissivity  factor  between  element  in  question  and 
surroundings.  is  evaluated  at  ^  (T^,T^). 

may  write  q^^  in  the  apparently  linear  form 


rs 


a  A  F  F  (T  ,T  )  P(T  ,T  )  (T  -T  ) 
A  ®  r  s  r  s  r  s 


where 


p  (T  T  )  =  T  +  T 

r  s  r  rs  rs  s 

In  this  way  we  may  obtain  from  the  formula 

K  =  <7  A  F.  F-  (T  ,T  )  P  (T  ,T  )  .  (9) 

rs  A  6  '  r  s  r  s 

Our  final  type  of  thermal  conductance  arises  from  forced  convec¬ 
tion.  Here  the  conductance  is  calculated  from  a  heat  transfer  coefficient 
h^g  and  is  given  by 


K 


rs 


A 


rs 


h 


rs 


(10) 


The  heat  transfer  coefficient  itself  is  calculated  in  a  more  complicated 
way  than  any  of  the  quantities  that  have  been  dealt  with  up  to  the  present. 
The  details  of  this  computation  are  described  in  Appendix  A  at  the  end  of 
this  report.  For  the  present  it  is  convenient  to  regard  it  as  a  function 
which  depends  on  time  and  on  the  tetrperature  T^.  In  the  application  of  a 
conductance  arising  from  aerodynamic  heating  in  the  heat  transfer  model  one 
end  point,  suppose  r,  of  each  such  conductance  has  the  applied  temperature 
T^aw.  The  calculation  of  this  temperature  is  described  in  Appendix  A  also. 


7. 
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In  applications  in  the  model,  the  heat  transfer  conductances 
usually  occur  along  with  conductances  arising  from  radiation  heat  transfer 
to  the  surrounding  space.  Here  the  radiation  transfer  has  the  special  form 


4  4 

q  -  a  A  6(T  )  (T  “)  -  T 

^rs  s  s  r  s 

With  the  conductance  given  in  the  corresponding  form: 


(11) 


K  «  A  a  6(T  )  P  (T  •  ,T  ) 
rs  s  s  r  s 


(12) 


here  is  a  temperature  that  corresponds  to  the  effective  temperature 

of  the  surrounding  space  and  e (T  )  is  given  in  tabular  form  or  by  a  poly- 

8 

nomial.  may  be  defined  as  a  function  of  altitude  or  time. 

III.  Specification  of  the  Model 

In  this  section  we  describe  in  some  detail  the  input  required  by 
the  computer  program  which  makes  it  possible  for  the  IBM  7090  computer  to 
perform  the  vast  amount  of  computation  involved  in  the  model  of  section  II. 

Other  questions  concerning  the  use  of  the  program  will  probably 
arise.  Some  of  these  questions,  such  as  those  relating  to  cost  (Machine 
Time)  and  the  size  of  problem  which  the  program  is  capable  of  handling 
(Machine  Storage) ,  can  be  answered  by  referring  to  Appendix  D.  Other  of 
these  questions  will  be  best  answered  by  establishing  direct  liaison  with 
BCC. 

Input  consists  of  a  simple  but  complete  FORTRAN  control  program 
prepared  by  the  user  and  a  set  of  FORTRAN-FAP  subroutines  (Nos.  11.02.01- 
11.02.12)  available  in  the  BCC  Library.  The  control  program  must  be  prepared 
on  standard  FORTRAN  coding  sheets  which  are  obtainable  from  BCC. 

There  is  a  considerable  amount  of  flexibility  in  the  method  of 
writing  the  control  program.  For  the  sake  of  simplicity  we  describe  a 
method  which  requires  a  minimal  knowledge  of  FORTRAN  and  leave  it  to 
experience  to  suggest  modifications. 
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The  structure  of  a  typical  control  program  is  as  follows: 
*  C0NTR0L  STATEMENTS 

DIMENS  I0N  STATEMENTS 

C  C0MMENTS 

LISTS 

CALL  SET 

1  CALL  TRAJ 

CALL  AMBATM 

CALL  F0RCER 

CALL  F0RCER 
CALL  F0R  ALT 

CALL  F0R  ALT 
CALL  AER0 


CALL  AER0 
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CALL  RAD 


CALL  RAO 
CALL  CdiM  C0N 

CALL  C0M  C^K 
CALL  C0N 

« 

CALL  C0N 

CALL  CAP 


CALL  CAP 

CALL  WRITE 

CALL  PL(^ 

CALL  PL^ 
CALL  STEP 

Gdi  T0  1 


END 
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The  effect  of  this  control  program  Is  to  supply  the  FORTRAN  MONITOR  and 
the  various  subroutines  with  the  Information  they  require  and  to  sequence 
the  machine  once  per  time  step  through  the  subroutines  In  the  order  In 
which  they  are  called  starting  with  statement  1. 

Reference  to  section  II  and  Appendix  A  should  clarify  the 
physical  meaning  of  all  the  parameters  that  occur  In  the  various  CALL  state¬ 
ments,  and  reference  to  section  V  should  clarify  the  form  In  which  they  are 
written  on  the  coding  sheets  as  well  as  the  Interplay  between  these  para¬ 
meters  and  the  LISTS.  The  definitions  of  all  parameters  In  the  various 
CALL  statements  follows: 

1.  CALL  SET  (START,  STOP,  TEMPIN,  I,  INDEXl,  Tl,  — ,  INDEXI  TI) 

START,  STOP  -  Start  computation  at  time  START  and  stop  computation  at 
time  STOP. 

TEMPIN, I, INDEXl, Tl,— INDEXI, TI  -  All  thermal  capacitors  are  Initialized 
to  TEMPIN,  with  exception  of  those  I  thermal  capacitors  whose  Indices 
are  listed.  INDEX  J  Is  the  Index  of  the  Jth  exceptional  thermal  capacitor 
which  Is  Initialized  to  temperature  TJ, 

2.  CALL  TRAJ  (FOFXM,XM,FOFXA,XA) 

NOTE:  The  pair  F0FX,X  occurs  In  CALL  TRAJ,  CALL  AMBATM,  CALL  F0R  ALT, 

CALL  F0RCER,  CALL  AER0,  CALL  RAD,  CALL  C0MC0N,  CALL  C0N,  and  CALL  CAP. 

It  Is  used  to  represent  a  function  FOFX(X)  In  one  of  the  following  two 
ways : 

(1)  to  define  FOFX(X)  as  the  polynomial 
N 

FOFX(X)  -  ^  AjX^ 

1-0 

write  N  for  X  In  the  CALL  statement  and  write  the  following  In  the 
data  LISTS: 

FOFX(l)  -  Aj^ 

F0FX(2)  - 

A 

o 


F0FX(1«-1) 
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(11)  to  define  FOFX(X)  as  a  table, write  the  following  In  the  data  LISTS: 
X(l)  -  N  . 

X(2)  -  X^ 


X(lW-l)  -  Xjj 
FOFX(l)  *  FOFX(Xj^) 


FOFX(N)  =«  FOFX(Xjj) 

The  values  FOFX(xy  are  .'theft  defined  by  linear '  Interpolation' 

In  CALL  AMBATM  and  by  quadratic  Interpolation  In  all  other  CALL 
statements. . 

FOFXM,XM  -  This  pair  represents  Mach  number  as  a  function  of  time. 
FOFXAjXA  -  This  pair  represents  altitude  as  a  function  of  time. 

3.  CALL  AMBATM  (FOFXTO,XTO,PO,N) 

FOFXTO,XTO  ^  This  pair  represents  ambient  temperature  as  a  function 
of  altitude. 

PO,N  -  Ambient  pressure  Is  determined  by  writing  In  the  data  LISTS. 
P0(1)  -  Zj  *  0 

P0(2)  -  ttj 

P0(3)  -  pj 

P0(4)  -  Yi 

P0(5)  - 

PO(6)  -  O2 

PO(7)  -  P2 

P0(8)  -  Y2 

PO(4N-3)  - 

PO(4N-2)  - 

PO(4N-l)  -  Pj^ 

P0(4N)  -  Yj, 


where  the  relationship  between  the  a,  p.Yand  Z  are  as  per  Equations  A5 
and  A6  of  Appendix  A. 
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4.  CALL  F0RCER  ( INDEX, FOFX.X) 

INDEX  -  This  is  the  Index  of  the  thermal  capacitor  whose  temperature 
Is  to  be  set  to  T. 

FOFX.X  -  This  pair  represents  T  as  a  function  of  time, 

5.  CALL  F0RALT  ( INDEX, FOFX.X) 

INDEX  -  This  is  the  index  of  the  thermal  capacitor  whose  temperature  is 
to  be  set  to  T. 

FOFX,X  -  This  pair  represents  T  as  a  function  of  altitude. 

6.  CALL  AER0  (ID,I,FL0N0S,F0FXLM,XLM,F0FXLP,XLP,INDEX1,GE0N01,P0SN01,  — , 
INDEXI , GEONOI , POSNOI) 

ID  -  This  is  aerodynamic  heat  transfer  block  identification  number. 

I  -  This  is  number  of  thermal  capacitors  in  this  block. 

FLONOS  -  This  is  list  of  5  numbers  which  describe  the  flow  for  this  block. 
They  should  be  presented  in  the  data  LISTS  as  follows : 

FLONOS (1)  =  Re 

cr 

FLONOS (2)  -  a 

FLONOS (3)  = 

FLONOS (4)  =  C^ 

FLONOS (5)  =  C 

“t 

FOFXUI,XU(  -  This  pair  represents  the  ratio  of  local  Mach  number  to 

ambient  Mach  number  as  a  function  of  ambient  Mach  number. 

FOFXLF,XLP  -  This  pair  represents  the  ratio  of  local  pressure  to  ambient 

pressure  as  a  function  of  ambient  Mach  number. 

INDEX J,GEONOJ,POSNOJ  -  This  triple  represents  the  Index  number,  surface 
area  A,  and  characteristic  length  X  respectively  of  the  Jth  thermal 
capacitor  In  this  block. 

This  CALL  AER0  format  must  be  supplied  for  each  heat  transfer  block. 

7.  GALL  RAD  (INDEXI, INDEX2,GE0N0,F0FX,X) 

INDEXI , INDEX2  -  This  pair  represents  the  indices  of  the  two  thermal 
capacitors  which  are  exchanging  heat  by  radiation. 

GEONO  -  This  Is  the  exchange  factor  multiplied  by  the  surface  area  A. 
FORX,X  -  This  pair  represents  emisslvlty, e,  as  a  function  of  temperature. 
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8.  CALL  C<»MC(8N  (INDEX1,INDEX2,GE0N01,GECN02,GE0N03,F0FX1,X1,F0FX2,X2, 
FOFX3,X3) 

INDEX1,INDEX2  -  This  pair  represents  the  Indices  of  the  two  thermal 
capacitors  of  this  composite  (contact)  conductance  pair. 
GE0N01,GE0N02,GE0N03  -  This  triple  represents  area  to  length  ratio 
(A/L)^,  area  A.  and  area  to  length  ratio  (A/L)^  respectively  proceeding 
from  the  thermal  capacitor  of  INDEXl  to  the  thermal  capacitor  of  INDEX2. 
FOFXl,Xl,FOFX2,X2,FOFX3  -  These  three  pairs  represent  respectively, 
the  conductivities  k  as  functions  of  temperature  proceeding  from  the 
thermal  capacitor  of  INDEXl  to  the  thermal  capacitor  of  1NDEX2. 

9.  CALL  CgN  (INDEXl, INDEX2,GE0N0,F0FX,X) 

INDEXl , INDEX2  -  This  pair  represents  the  Indices  of  the  thermal  capacitors 
of  this  conductance  pair. 

GEONO  -  This  represents  the  area  to  length  ratio  A/L. 

F0FX,X  <•  This  pair  represents  conductivity  k  as  a  function  of  temperature. 

10.  CALL  CAP  (INDEX, GEONO, FOFX,X) 

INDEX  •  This  Is  the  Index  number  of  this  thermal  capacitor. 

GEONO  -  This  Is  volume  V. 

F0FX,X  -  This  pair  represents  specific  heat  pc^  as  a  function  of 
temperature. 

11.  CALL  WRITE  (I,TIME1,TIME2,  — ,TIMEN) 

I  -  This  Is  the  total  number  of  complete  printouts.  A  complete  print¬ 
out  consists  of  all  aerodynamic  and  trajectory  Information  and  the 
temperatures  on  all  thermal  capacitors  and  contact  surfaces  which 
differ  from  TEMPIN.  Such  printouts  are  meant  for  gross  checks  and 
should  be  used  sparingly. 

TIME1,TIME2,  —  ,TIMEN)  -  The  nuinbers  1,2,-— ,N  represent  the  times  at 
which  complete  printout  will  occur. 

12.  CALL  PL0T  (ID,DELTEM,I,INDEX1, - , INDEXl) 

ID  -  This  Is  the  Identification  number  of  this  PL0T  group. 

DELTEM  -  A  card  Is  punched  each  time  the  temperature  of  the  thermal 
capacitor  whose  Index  Is  INDEXl  changes  by  DELTEM.  The  card  deck  result¬ 
ing  from  a  complete  run  may  be  sorted  on  ID  (columns  1-4)  and  then  tabu¬ 
lated  to  obtain  a  printout  of  X-Y  plotted  to  obtain  a  graph  of  temperature 
(columns  9-12,  13-16,-— ,69-72)  versus  time  (columns  5-8). 
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I  -  This  is  the  number  of  thennsl  cspacltors  in  this  PL0T  group. 

INDEXl,  —  jINDEXI  -  'The«eIljtdlce8:.'rcpT*««tit  the  Indices  of  the 
thermal  capacitors  whose  temperatures  are  to  be  punched. 

13.  CALL  STEP 

This  Is  a  sequence  Instruction  for  the  computer  to  re-start,  timewise, 
with  the  solution  of  the  next  set  of  finite  difference  equations. 

IV.  Printout  Format 

The  sheets  of  printout  produced  In  the  problem  run  contain  the 
history  and  principal  results  of  the  analysis.  Here  are  recorded  the 
capacitance  tenq)eratures  designated  for  printout  and  the  contact  tempera¬ 
tures  that  appear  In  those  composite  resistances  that  hhve  been  designated 
for  printout. 

Whenever  aerodynamic  heating  Is  Included  In  a  heat  transfer  model 
It  Is  essential  to  have  some  record  of  the  principal  quantities  which  in¬ 
fluence  this  mode  of  heat  transfer.  To  provide  this  Information  the  values 
of  the  following  quantities  are  provided  at  each  printout  time:  missile 
altitude  z,  missile  Mach  number  M^,  ambient  pressure  p^,  and  anibient  tenq>era- 
ture  T^.  In  addition  to  this,  characteristic  Information  Is  printed  out  for 
each  heat  transfer  block  (CALL  AER0  statement) ;  here  the  recorded  quantities 
are  the  local  Mach  number  M,  local  pressure  p,  local  air  teiig>erature  T, 
local  Reynolds  number,  adiabatic  wall  temperature  T  ,  and  the  heat  transfer 

•  clW 

coefficient,  normalized  by  replacing  x,  the  characteristic  length,  by  1;  l.e. 

1-a 

to  get  h,  divide  the  values  printed  but  by  x  . 

The  program  also  has  a  subroutine  available  which  enables  output 
capacitance  temperatures  to  be  plotted  as  functions  of  time  on  BID  X-Y 
plotters.  The  complexity  of  a  given  thermal  problem,  operator  experience, 
time  available  and  computer  budget  money  available  dictate  the  best  use 
that  can  be  made  of  the  program.  To  help  clarify  the  discussion  a  8aiq>le 
problem  Is  Included  In  the  next  section. 
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V.  EXAMPLE 


Figure  1  shows  a  wedge  type  airfoil  made  of  stainless  steel 
It  is  seen  that  the  airfoil  has  three  webs.  These  webs  are  made  of  stain¬ 
less  steel  #347.  They  are  of  0.051  inch  thickness.  The  skin  thickness  of 
the  airfoil  Is  0.040  Inch.  It  is  assumed  in  this  problem  that  the  vertex 
angle  of  the  wedge  is  4  degrees  and  that  one  surface  is  at  a  flxec  angle 
of  attack  of  4  degrees  with  the  alrstream,  while  the  other  surface  is  in 
line  with  the  alrstream.  Temperature  distributions  are  desired  on  all  skin 
surfaces  and  webs  during  a  prescribed  flight  pattern.  It  is  recommended 
that  the  reader  uses  the  FORTRAN  sheets  at  the  end  of  this  report  in  order 
to  follow  through  this  sample  problem. 

Geometry 

Figures  2  and  3  show  a  convenient  method  of  dlsecting  the  physical 
airfoil  into  a  lumped  parameter  system  of  capacitors. 

TRAJECTORY  FUNCTIONS 

To  define  the  trajectory,  the  altitude  and  free  stream  Mach 
number  are  prescribed  as  functions  of  time  under  the  TRAJECTORY  FUNCTIONS 
heading,  page  2.  They  appear  as  functions  FI  and  F2. 

AMBIENT  ATMOSPHERE  FUNCTIONS 

It  is  seen  in  Appendix  A,  Eq.  (A3) ,  that  in  order  to  evaluate  the 
expression  h,  the  local  flow  conditions  (pressures  and  velocities  that  exist 
immediately  outside  of  the  boundary  layer),  must  be  determined.  These,  in 
turn  depend  on  the  flight  trajectory  of  the  airfoil  and  on  the  distribution 
of  ambient  pressure  and  tenqjerature  with  altitude.  These  last  quantities 
are  referred  to  under  the  AMBIENT  ATM0SFHERE  FUNCTI0NS  heading.  Function  F3 
describes  the  variation  of  ambient  temperature  with  altitude.  The  function 
F4  describes  constants  to  be  used  in  the  computation  of  ambient  pressure  by 
functions  of  the  form  given  by  Eq.  (A5)  of  Appendix  A.  The  order  of  listing 
these  constants  is  as  follows. 


1st , 
2nd, 
3rd, 
4th, 


list  lowest  altitude  of  the  trajectory 
■  ^ 

list  a 

, which  satisfy  Eq.  (A5)  for  ambient  pressure  in  the  lowest 
to  intermediate  altitude  range. 


f 
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5th,  list  an  Intermediate  altitude  of  the  trajectory 
6th,  llst^ 

7th  list  p>  satisfy  Eq,  (A5)  for  ambient  pressure  in  the  lntermedlat( 


8th,  list  Y 


to  high  altitude  range. 


9th,  list  a  high  altitude  of  the  trajectory 
10th,  list  gl 

11th  list  satisfy  Eq,  (A5)  for  ambient  pressure  above  this  high 


12th,  list  Y 


altitude 


FLOW  NUMBERS 

As  may  be  seen  In  Eq.  (A3) ,  In  order  to  evaluate  the  coefficient  of 
convective  heat  transfer  h,  the  constant  Cq  and  the  exponent  0l  have  to  be 
specified.  The  constant  is  dependent  on  the  type  of  flow  existing  on  the 
surface  of  Investigation,  l.e. ,  flat  plate,  duct  or  cone,  and  on  the  transi¬ 
tion  Reynolds  number  which  Indicates  when  the  flow  changes  from  laminar  to 
turbulent.  The  exponent  Is  solely  dependent  on  the  transition  Reynolds  number. 
Consequently,  under  the  FL0W  NUMBERS  heading,  the  transition  Reynolds  number 
Is  specified.  It  Is  followed  in  order,  by  the  value  of  a  for  Reynolds  number 
less  than  or  equal  to  the  critical  value,  the  value  of  a  for  Reynolds  number 
greater  than  the  critical  value,  the  value  Cq  for  Reynolds  number  less  or 
equal  to  critical,  and  the  value  of  C^  for  Reynolds  number  greater  than  critical. 
In  our  sample  problem  we  consider  flat  plate  theory  in  which  we  assume  R^~1.5xl0^, 
a  ■  0.5  laminar,  o«0.8  turbulent,  Cq  -  0.332  laminar,  and  Cq  -  0.0265  turbulent. 
The  above  constants  and  exponents  appear  as  functions  F5. 

L0CAL  FL0W 

The  local  flow  conditions  are  obtained  by  the  presentation  of  the 
ratio  of  local  pressure,  p,  to  the  free  stream  pressure,  p^,  and  ratio  of 
local  Mach  number,  M,  to  the  free  stream  Mach  number,  M^,  as  functions  of 
free  stream  Mach  number,  M^.  In  the  case  of  elements  1,  2,  3,  4,  5  of  the 
airfoil  It  is  assumed  that  the  local  Mach  number  and  the  local  pressure  are 
the  same  as  the  respective  free  stream  Mach  number  and  free  stream  pressure. 
Therefore,  the  ratio  of  unity  appears  as  function  F6.  For  the  other  capaci¬ 
tances  the  ratios  are  presented  In  the  form  of  polynomial  functions  F7  and 

F8.  The  form  M  (or  p)*f(M)»ax’*  +  a  ,x*'~^  ••••  a„x^  +  a.x  +  a  Is 

'o'  n  n-1  2  1  o 
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used  in  which  the  coefficients  a  ,a  ***  a  are  such  that  the  continuous 

o  1  n 

function  f(M  ),  closely  circumscribes  the  discrete  points  M/M  (and  p/p  ) 
O  0  0 

for  which  experimental  data  Is  available.  The  ratios  M/M^  and  p/p^  can  be 

listed  In  tabular. form  as  functions  of  M  If  sufficient  data  points  are 

o 

available.  In  that  case,  a  quadratic  Interpolation  Is  used  between  data 
points. 


EMISSIVITIES 

Under  this  heading,  the  effective  emlsslvlty  of  the  exterior  skin 
of  the  airfoil  Is  listed  as  F9.  In  the  sample  problem  It  Is  assumed  that 
the  emlsslvlty  remains  at  the  constant  value  of  0.8  at  all  skin  temperatures. 

THERMAL  C0NDUCTIVITY  and  THERMAL  HEAT  CAPACITANCE 


As  was  Indicated  In  Section  II,  the  thermal  properties  of  each 
capacitance  and  conductance  must  be  specified  as  functions  of  tenqperature. 
For  the  properties  of  this  sample  problem,  the  thermal  conductivity  k,  and 
the  thermal  heat  capacitance  pC^  are  listed  as  polynomial  functions  of 
temperature  of  1st  and  3rd  degree,  respectively.  They  appear  as  functions 
FIO  and  F12,  respectively,  on  page  6  of  the  FORTRAN  sheets.  These  poly¬ 
nomials  were  selected  by  examining  the  proximity  of  the  functions  with  the 
discrete  points  from  which  they  were  generated.  In  the  past,  these  poly¬ 
nomials  were  generated  on  the  IBM  650  computer.  At  present  BBE  has  a 
library  of  such  polynomials  for  a  great  variety  of  metals  and  non-metals. 

CALlI  SET  Statement  (F0RTRAN  page  7) 


In  order  for  the  computations  to  be  carried  out  It  Is  necessary 
to  specify  the  flight  time  at  which  to  Initiate  the  program  and  the  flight 
time  at  which  to  cease  such.  This  appears  as  0  and  100,  respectively  In 
the  CALL  SET  statement.  This  Information  Is  followed  by  the  Initial 
temperature  that  Is  assigned  to  each  capacitance.  In  the  present  case 
this  appears  as  560.  If  there  is  radiation  by  external  surfaces  to  space, 
as  in  the  present  case.  It  Is  necessary  to  list  the  temperature  that  Is  to 
be  assigned  to  space.  In  the  present  case  capacitor  1000  refers  to  space, 
and  Its  temperature  appears  as  517.  If  the  temperature  of  space  varies  with 
altitude  (or  flight  time).  It  Is  necessary  to  state  the  Initial  temperature 
of  space  In  the  CALL  SET  statement  and  In  addition.  It  Is  necessary  to  list 
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the  variation  of  the  space  temperature  either  with  altitude  or  with  time. 

The  former  would  be  called  out  by  the  CALL  F0RALT  subroutine,  while  the 
latter  by  the  CALL  F0RCER  subroutine.  See  parts  4  and  5  of  Section  III. 

CALL  TRAJ  Statement 

The  CALL  SET  statement  is  followed  by  the  CALL  TRAJ  statement. 

The  latter  calls  out  the  Mach  No.  FI  vs.  time  FIT,  and  altitude  F2  vs. 
time  F2T  functions. 

CALL  AMBATM  Statement 

The  CALL  AMBATM  statement  follows  the  CALL  TRAJ  statement.  It 
calls  out  the  ambient  temperature  F3  vs.  altitude  F3Z  functions  and  the 
polynomial  function  F4  which  describes  the  pressure  as  a  function  of 
altitude.  The  degree  of  the  polynomial  is  listed  after  the  polynomial 
function  F4  in  the  CALL  AMBATM  statement.  This  is  necessary  in  order  for 
the  computer  subroutine  to  introduce  all  of  the  polynomial  coefficients  into 
the  program  run. 

CALL  AER0  Statement 

The  CALL  AER0  statements  specify  the  computations  for  the  con¬ 
ductances  and  adiabatic  wall  temperatures  which  arise  from  aerodynamic  heat¬ 
ing.  Each  statement  describes  a  heat  transfer*  block,  consisting  of  separate 
conductances  driven  by  the  same  adiabatic  wall  temperature,  computed  from 
the  same  local  flow  conditions  and  reference  temperature  T' ,  radiating  to 
the  same  space  temperature,  but  having  different  areas,  reference  lengths, 
and  emlssivlties.  Each  block  is  given  a  block  number.  Blocks  are  numbered 
in  consecutive  order.  Adjacent  to  the  block  number  In  the  CALL  AER0:  state¬ 
ment  is  a  number  which  states  how  many  conductances  are  described  by  the 
block.  Next  to  this  number  the  function  F5  from  the  FL0W  NUMBER  heading 
appears.  This  is  followed  by  the  pertinent  L0CAL  FL0W  functions  F6  or  F8  and 
F7.  It  is  Important  to  note  that  the  function  describing  the  local  Mach  No. 
has  to  be  listed  first,  and  is  always  followed  by  that  function  which  describes 
the  local  pressure.  Again,  these  functions  are  followed  by  the  degree  of  the 
polynomial.  Next,  the  most  representative  element  -  element  6,  in  the  sense 
of  thickness  and  material  properties  is  listed.  This  is  followed  by  the 
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aerodynamic  surface  area  -  0.0416  ft  ,  and  in  turn  by  the  characteristic 
length  -  0.1408  ft.  Then  elements  1,  2,  3  etc.  are  defined  until  all 
elements  of  the  heat  transfer  block  are  described. 

The  temperature  of  that  capacitance  which  appears  first  in  the 
CALL  AER0  statement,  is  used  as  the  Tw  in  Equation  A2  to  evaluate  T'  for 
the  entire  block.  Otherwise,  the  order  in  which  the  capacitances  are 
listed  within  the  block  is  of  no  importance. 

CALL  RAD 

A  sample  page  of  the  CALL  RAD  statement  is  shown  on  F0RTRAN 
Sheet  8.  Here,  the  radiation  interchange  between  exterior  skin  surfaces 
and  space  is  considered.  First  is  listed  the  element,  then  the  identifica¬ 
tion  number  for  space,  then  the  element  surface  area.  In  turn,  the  emisslvlty 
function  F9  is  listed,  and  in  turn,  is  followed  by  the  degree  of  this  poly¬ 
nomial  which  is  zero. 

CALL  GtSACm 

A  sample  page  of  the  CALL  C0MC0N  Statement  is  shown  on  F0RTRAN 
Sheet  9.  In  the  region  where  heat  conduction  from  the  surface  elements 
to  the  webs  is  considered,  the  composite  conductance  technique  described 
in  Section  II  is  applied.  The  identifying  elemental  numbers  are  listed  as 
1,2,  .  Then,  the  A/L  ratios  are  listed  in  the  following  order:  first, 

9 

that  A/L  ratio  which  pertains  to  element  1,  then  a  large  number  1  x  10  to 
make  the  contact  resistance  0,  then  the  A/L  ratio  which  pertains  to  element 
2.  In  turn,  the  polynomial  function  FIO  which  describes  the  thermal  con¬ 
ductivity  of  the  material  of  element  1  and  its  degree  are  listed.  This  is 
followed  by  the  number  1  which  is  a  multiplier  for  the  large  number  previously 
mentioned.  Thus,  the  product  is  a  large  number  and  this  makes  the  contact 
resistance  0.  Lastly,  the  polynomial  function  describing  the  thermal  con¬ 
ductivity  of  element  2,  and  its  degree  are,  respectively,  listed. 

CALL  C0N 

An  illustrative  sheet  of  the  CALL  C0N'  statements  is  shown  on 
F0RTRAN  Sheet  10.  The  order  of  listing  these  statements  should  now  be 
apparent.  The  two  identifying  capacitance  numbers  are  listed,  and  are 
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followed  by  the  A/L  tern,  and  In  turn,  by  the  function  describing  the 
thermal  conductivity,  and  the  degree  of  the  polynomial. 

CALL  CAP 

An  Illustrative  sheet  of  the  CALL  CAP  statements  Is  shown 
F0RTRAN  Sheet  11.  First  Is  listed  the  capacitance  number,  then  the 

3 

volume  In  cu.ft. ,  then  the  heat  capacitance  pC^  (In  the  units  of  BTU/ft  -“R) 
F12  polynomial  function,  and  finally  the  degree  of  the  polynomial. 

CALL  WRITE 

Here  are  written  the  flight  times  at  which  the  temperatures  of 
the  capacitors  Is  to  be  printed  out.  First,  Is  listed  the  total  number 
of  time  steps  33  for  which  prlnteout  Is  requested,  then  the  times  at 
which  prlnt-out  Is  desired  are  stated,  l.e.,  0.,  1.5,  etc. 

CALL  PL0T 

A  subroutine  exists  which  enables  output  temperatures  to  be 
punched  on  cards  which.  In  turn,  may  be  used  to  obtain  X-Y  plots.  See 
part  12,  Section  III.  This  procedure  Is  not  described  here  as  It  has 
been  found  more  convenient  to  plot  the  temperatures  by  hand. 

CALL  STEP.  G0  to  1.  END 

This  Is  the  Information  necessary  for  the  computer  to  continue 
timewise  with  the  solution  of  the  next  set  of  finite  difference  equations, 
and  continues  until  the  expiration  of  the  flight  time. 

Results 

Figure  4  shows  a  typical  prlnt-out  page  of  this  sanple  problem 
with  the  proper  Identifications.  Though  these  temperatures  appear  In 
degrees  Ranklne,  the  program  has  recently  been  modified  so  that  tempera¬ 
tures  now  are  printed  out  In  degrees  Fahrenheit.  The  pressure  Is  expressed 
In  pounds  per  square  foot,  and  the  altitude  Is  expressed  In  feet.  The 


TH«  iohni  Hopliini  Unlvtrtity 

AFPtlSO  PMVtiet  UMRATPRV 

SilvAr  Spring,  MUrylAnd 


21 


1  “(X 

normalized  heat  transfer  coefficient  h'  when  divided  by  x  is  expressed 
2 

in  BTU/ft  -sec-^R  units. (This  is  the  convective  heat  transfer  coefficient  h.) 

The  results  of  this  problem  are  presented  in  the  form  of  plots 
of  temperature  as  a  function  of  distance  from  the  leading  edge  at  various 
Instants  of  time  for  both  sides  of  the  wedge.  In  addition,  a  plot  of 
temperature  vs.  time  is  presented  across  one  of  the  webs.  These  plots 
are  shown  in  Figures  5,  6  and  7.  It  may  be  seen  that  temperature  drops 
occur  for  the  surface  elements  close  to  the  web.  This  is  due  to  the  heat 
sink  effect  of  the  webs. 

It  is  hoped  that  the  above  sample  problem  has  sufficed  in  intro¬ 
ducing  the  reader  to  the  BBE-BCC  heat  transfer  program.  A  more  comprehen¬ 
sive  understanding  will  come  about  in  the  actual  usage  of  the  program. 
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APPENDIX  A 

In  the  main  body  of  this  memorandum  we  have  referred  to  a 
calculated  heat  transfer  coefficient  and  accompanying  driving  tempera¬ 
ture.  This  appendix  sets  forth  the  manner  In  vdilch  these  quantities 
are  calculated. 

The  heat  transfer  coefficients  that  are  used  are  based  on 
a  formula  of  the  form 


®  P  k 

h  -  Ctt  Re  5  .  (Al) 

In  which  la  a  constant,  Is  a  Reynolds  number,  P^  is  the  Prandtl 
number  for  air,  k  Is  the  thermal  conductivity  of  air,  and  x  Is  a  length. 

In  our  applications  we  will  compute  these  quantities  from  "local  flow" 
pressures  and  velocities,  l.e. ,  those  which  exist  Immediately  outside  the 
boundary  layer,  and  from  an  effective  temperature  T'  given  by  the  formula 


T'  •  Ma  \  +  ''o  (1  +  ®o  T 


(A2) 


where  Pg, v^,  and  Oa  are  constants,  M  and  T  are  local  Mach  number  and 
teiq>erature  respectively,  and  T^  is  the  temperature  of  the  surface  to 
which  aerodynamic  heat  transfer  takes  place. 


By  using  the  perfect  gas  law  and  associated  specific  heat 
relationships 


C 

P 


£R 

Y-1  J  ’ 


the  definition  of  Reynolds  number  and  Prandtl  number 
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the  Sutherland  expression  for  the  viscosity  of  air  as  a  function  of 
temperature. 


and  the  expression  (A2) ,  the  values  of  h  can  bs  obtained  from  expression 
(Al)  as  follows, 


and 


R 

e 


P  vx 


p  V  x 
gRT'  M> 


k 


P  P  Y-1  J 
r  r 


Therefore, 


h 


4.  Cl  Op  1  .  -  «Ct  >  , »  “O  1  “Ct 

(pv)  X  ^(gR)  (T')  \i 


18-1 

,(T’) 

. 


-1  £R 
J 


Y  (T’) 
Y(T')-1 


By  collecting  terms. 


h 


(gR)^‘®(pv)V’V^‘“(T’)' 


1. 


(T') 


P-1 


Y  (T') 

y(t')-i 


Now,  by  Introducing  the  Sutherland  expression  the  expression  for  h 
becomes , 


-a 


, 

(T.)3/2 

1-a 

(T’)"“j’p^(T’) 

, 

T'+T 

c 

L 

" 

P-1 


Y  (T') 


y(t')- 


(A3) 


Here  C^j,  J,  |a.^,  g  and  R  are  constants,  p  Is  the  local  pressure,  v  Is 
the  local  velocity,  x  Is  a  length,  T'  Is  given  by  (A2) ,  Is  a 
constant  temperature,  and  Y Is  the  ratio  of  specific  heats  for  air.  In 
the  expression  (A3)  It  may  be  observed  that  the  entire  quantity  within 
.the  braces  Is  a  function  of  T'  alone  and  Its  value  Is  completely  determined 
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when  the  variations  of  and  Y  functions  of  temperature  are  prescribed. 

Here  we  will  calculate  y  ^nd  from  approximations  by  fifth  degree  poly¬ 
nomials  over  the  range  of  interest. 

In  order  to  evaluate  the  expression  for  h,  the  local  flow  condi¬ 
tions  must  be  determined.  These  in  turn  depend  on  the  trajectory  of  the 
missile  flight  and  on  the  distribution  of  ambient  pressure  and  temperature 
with  altitude. 

We  will  describe  the  variation  of  ambient  temperature,  T^,  with 
altitude,  z,  as  a  continuous  piecewise  linear  function  of  altitude,  i.e.. 

To  -  5^  +  z,  z^_j  <  z  <  z^,  (i»l,2,‘** ,5) ,  (A4) 

where  <>,  aad  k  are  constants.  Distributions  of  ambient  pressure,  P  , 

11  o 

hydrostatically  compatible  with  linear  temperature  distributions  are  of 
one  of  two  forms ,  either 

Y. 

Po  "  0^11+  (A5) 

i  or 

Pq  =  exp  [P^  (z  -  z^_j)  ]  ,  zj^  j  <  z  <  z|  .  (A6) 

The  first  is  compatible  with  a  non-constant  linear  temperature  distribution 
from  z|  to  zj^,  while  the  second  corresponds  to  an  isothermal  temperature 
distribution  over  the  same  interval. 

To  define  a  trajectory  with  sufficient  accuracy  for  most  heat  trans¬ 
fer  purposes  the  altitude  z  and  free  stream  Mach  number,  M^,  must  be  prescribed 
as  functions  of  time.  These  are  given  as  continuous  piecewise  linear  func¬ 
tions  of  time.  Such  a  function  is  completely  defined  when  the  coordinates 
are  given  for  its  initial  point,  for  its  final  point,  and  for  all  "corners" 
at  which  the  slope  of  the  function  changes. 
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The  local  flow  conditions  can  then  be  obtained  provided  that  the 
ratio  of  the  local  Mach  number,  M,  to  the  free  stream  Mach  number,  M  and 
the  ratio  of  (be  local  pressure  p,  to  the  free  stream  pressure,  p^,  are 
prescribed.  Here  we  will  list  these  In  tabular  form  as  functions  of  free 
stream  Mach  number. or  as  polynomial  functions  of  M^. 

The  remainder  of  the  local  flow  conditions,  l.e.,  the  local 
temperature,  enthalpy,  and  velocity  are  to  be  computed  from  the  relations 

H(T)  +  It  T  Y(T)  M^  -  H(T  )  +  |t  T  Ytt)  M  ^  (A7) 

/J  O  o  o  o 


and 


V 


[g  RT  Y(T)  M^ 


(A8) 


Here  H  Is  the  enthalpy  and  R  the  "gas  constant"  for  a  pound  of  air  and  J  Is 
the  conversion  factor  from  energy  In  mechanical  units  to  energy  In  thermal 
units.  We  will  approximate  H(T)  over  the  range  of  Interest  by  a  polynomial 
of  fifth  degree. 

Returning  to  the  expressions  for  the  heat  transfer  coefficient  and 
T' ,  It  will  be  noticed  that  several  constants  have  been  given  the  subscript 
a,  the  exponent  In  (Al) .  The  purpose  of  this  notation  Is  to  Indicate  that 
during  the  course  of  the  confutations  a  may  be  changed,  and  along  with  It 
the  constants  indexed  by  It.  In  particular  we  will  change  a  when  the  local 
Reynolds  number  exceeds  a  specified  critical  value  associated  with  transition 
from  laminar  to  turbulent  flow.  The  local  Reynolds  number  Is  calculated 
from  the  formulas 


R 

e 


pv  X 


(A9) 


P 


_2 _ 

R  g  T 


t3/2 

**0  T  +  T  ’ 
c 


(AlO) 

(All) 


where  p,  v,  and  T  are  local  flow  quantities. 
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The  temperature  which  acts  as  the  source  or  driving  temperature 
for  aerodynamic  heat  transfer,  l.e. ,  the  "adiabatic  wall"  temperature,  is 
computed  on  an  enthalpy  basis  as  follows: 


-  H  (T)  +  r^  T  Y  (T)  m2 


T  -  T  (H  ) 
aw  aw. 


(A12) 


in  which  r^  is  a  recovery  factor  and  T(H)  is  the  function  which  gives  air 
temperature  as  a  function  of  enthalpy.  This  last  function  can  be  approxi¬ 
mated  over  the  range  of  interest  by  a  polynomial  of  degree  five. 
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APPENDIX  B 


Stability 

While  there  Is  available  no  exact  treatment  of  stability  for  the 
type  of  equations  that  we  consider  In  our  heat  transfer  model,  It  Is  reason¬ 
able  to  expect  (and  experimental  results  bear  this  out)  that  the  same  criterion 
that  Is  sufficient  for  stability  In  those  cases  which  can  be  treated  exactly 
will  be  suitable  here.  Briefly,  this  criterion  Is  that  the  Inequality, 


At 


(Bl) 


hold  for  all  points  j  which  have  capacitance  associated  with  them.  Here, 
of  course,  the  sum  Is  taken  over  all  Indices  1  which  are  connected  to  j 
by  non-zero  capacitances. 


This  Inequality  may  be  given  the  following  physical  Interpretation: 
If  the  values  of  the  conductances  and  capacitances  In  the  network  are  calcu¬ 
lated  for  any  point  j  for  the  actual  temperature  distribution  that  occurs  In 
the  model  at  a  given  time,  then  the  time  step  allowable  at  that  time  Is  such 
that  even  If  the  surrounding  tenqjeratures  were  set  equal  to  zero  the  heat  that 
would  flow  out  of  Cj  In  the  time  At,  l.e.  ,  ^  Tj  At,  would  not  be 

sufficient  to  cause  the  temperature  Tj  +  A’^j  to  become  negative.  In  this 
form  the  criterion  (Bl)  may  be  regirdtd  as  an  expression  of  the  second  law 
of  thermodynamics.  In  f  ur  case  the  quantities  and  are  functions  of 
temperature  and  are  thus  not  known  before  the  analysis  Is  completed.  This 
means  that  In  order  to  evaluate  the  stability  of  the  problem,  we  must  make 
some  a  priori  estimates  of  the  tesqierature  distributions  that  are  expected. 


The  machine  program  has  been  set  up  so  that  If  we  set  START  time  = 
STOP  time  In  the  CALL  SET  statement  a  one  loop  computation  will  result  and 
a  complete  printout  plus  two  lines : 


CRITICAL  INDEX  Np.  »=  I  III 
DELTIM  -  0.  XXXXXXXXE  ±  XX 

will  be  Indicated.  A  hand  conq>utatlon  of  At  for  thermal  capacitor  IIII 
should  now  be  made  with  equation  (Bl).  If  there  Is  more  than  a  107.  difference 
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between  this  hand  calculated  A  t  and  the  DELTIM  computed  by  the  machine  a 
non-dlagnostlc  error  In  the  control  program  exists. 

An  estimate  of  the  number  of  time  steps  required  to  compute  from 
START  to  STOP  Is  furnished  by 

_  STOP-START 
DELTIM 
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APPENDIX  C 

BBE-BCC  Finite  Difference  Heat  Transfer  Equations 


The  Fourier  differential  equation  which  expresses  the  rate  of 
flow,  or  flux,  of  heat  through  a  given  surface  element  In  the  direction 
of  Its  normal  Is  proportional  to  the  normal  derivative  of  the  temperature; 
that  Is, 


<P  - 


(1) 


where  dA  Is  the  area  of  the  surface,  and  k  Is  the  conductivity.  Then  If 
V  Is  a  volume  enclosed  by  a  surface  S,  the  rate  of  gain  of  heat  by  the 
volume  V  is 


<P  - 


kVTdV 


where  the  operator^^  (read  "del") 


Is  defined  as 


V- 


(2) 


and 


is  a  scalar  product  called  the  divergence  of  1\/t. 


But,  If  p  Is  the  density  and  the  specific  heat  of  the  material  Involved, 
the  rate  of  gain  of  heat  Is  also  given  by 


and  so 


9 


d  V  , 


(3) 


9  (pCp)T 
3t 


d  V 


0  . 


(4) 


V 
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But  this  Is  true  for  any  volume  V,  so  the  integrand  Itself  must  vanish  for 
every  point  In  the  material;  thus  we  arrive  at  the  heat  equation, 

d(pC  )T 

^^7*^  “at 


This  Is  a  partial  differential  equation  which  governs  the  flow  of  heat  by 
conduction  In  a  homogeneous  or  non-homogeneous ,  Isotropic  or  anisotropic 
medium  which  Is  free  of  sources  and  sinks.  One  method  of  approximating 
Equations  (5)  or  (6)  for  conduction  In  one-dlmenslon  Is  by  the  following 
finite  difference  equation. 
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Here  for  a  finite  time  Interval  At,  the  term  AT/Ax  Is  expressed  in  terms  of 
conditions  at  the  beginning  of  the  time  interval.  Therefore,  this  difference 
equation  is  called  a  forward  difference  equation.  (Central  and  backward 
difference  equations  are  also  possible  to  construct,  but  will  not  be  the 
subject  of  this  Appendix.) 

It  is  seen  that  with  the  above  equation  the  temperature  at  any 
point  X,  can  bp  found  at  a  future  time  t+At  provided  only  that  a  history  of 
the  present  temperature  of  the  point  x,and  neighboring  points  x^Ax  and 
x-Ax  is  known.  However,  this  equation  holds  only  for  conduction  in  a 
medium  which  is  free  of  sources  and  sinks.  If  the  desired  internal  tempera¬ 
ture  is  that  of  a  point  which  is  adjacent  to  a  boundary  point,  i.e. ,  point 
x^  next  to  boundary  point  x^,  then  it  is  necessary  to  know  the  temperature 
of  point  x^  at  a  previous  time.  For  example,  in  order  to  evaluate  the 
temperature  of  point  X2  at  time  t^  by  Equation  (7),  it  is  necessary  to  know 
the  temperature  of  element  Xj^  at  time  t^.  This  is  found  by  rewriting 
Equation  (1)  in  the  following  form. 


conv' 

_  cont.res. 

where  now  refers  to  the  heat  applied  at  the  surface  of  Xj^  by  convection 

and/or  radiation  or  contact  with  some  other  medium  i.e.,  contact  resistance. 
Writing  Equation  (8)  in  a  finite  difference  form  leads  to  the  following, 


E 


rad. 

conv. 

cont.res. 


q  (xj.t^.t^)  *  -  k  (T^  Xj) 


^^x  t  ■  ’^x  t  ^ 
^2*  2  1*  2 

A  X 


(9) 


If  y  q  is  known  then  the  only  unknown  parameter  in  Equation  (9)  is  T  . 
Therefore,  T  ^  evaluated.  It  is  then  used  in  Equation  (7)  to  evaluate 
T  at  t«.  If  Vq  is  not  known,  it  is  always  expressed  in  a  form  in  which 
it  is  linearly  dependent  on  T  .  Therefore,  Equation  (9)  can  still  be 

Xi,t2 

solved  for  T  .  This  linearization  is  somewhat  demonstrated  by  Equations 
(5),  (6),  (9),  (10),  and  (12)  of  the  text.  The  additional  necessary  Information 
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is  that  in  the  case  of  a  non-linear  equation  such  as  Equation  (8)  of  the 
text,  the  non-linear  part  K  ,  is  always  evaluated  at  a  previous  time  Instant, 

ITS 

i.e.,  in  the  above  discussion  K  would  be  evaluated  at  t,. 

*  rs  1 
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Dimensional  Units 


Since  the  quantities  involved  In  the  heat  transfer  calculations 
Involve  dimensional  quantities  It  Is  necessary  to  establish  the  dimensional 
units.  Throughout,  the  English  system  of  units  Is  used,  l.e. ,  foot,  pound, 

°R,  BTU,  etc.  Time  by  Itself  Is  always  In  seconds,  but  when  stating  thermal 
quantities  such  as  conductivity  we  shall,  for  convenience,  use  hours. 

The  following  table  lists  the  units  used  for  all  quantities  that 
must  be  listed  on  the  problem  sheets  or  appear  on  the  printout  sheets.  (Since 
the  completion  of  the  Included  sample  problem  a  modification  has  been  made  In 
the  computer  program  whereby  temperatures  are  printed  out  In  degrees  Fahrenheit 
rather  than  degrees  Rankine.  However,  Input  temperatures  still  have  to  be  listed 
in  *R.) 

Table  of  Dimensional  Units 


Quantity 
Initial  time,  t 

o 

Initial  Values 
Time  step,  t 
Atmosphere: 


'*1 

Pi 

^1 

Z. 


Specific  Heat,  pc^ 

Volume ,  V 
Area ,  A 

Area/Length,  A/L 
Conductivity,  k 
Stefan  Boltzmann  constant.  O’ 
Exponents ,  a 

Heat  Transfer  Constants,  C^ 
Characteristic  Length,  x 
Heat  transfer  Coefficient,  h 
Temperature,  T 
Geometric  exchange  factor, 

Net  emlssivlty  factor,  Fg 


Dlmens Ion 
second 
‘R 

second 


lb  ft"^ 


non-dimensional 


ft 

“R 

'R  ft 


ft 

BTU  ft' 
ft^ 


ft‘ 


ft 

BTU  hr"^  ft"^  “r"^ 
BTU  hr"^  ft"'^  “R"^ 


non-dlmens lonal 
non-dimensional 
ft 

BTU  sec"^  ft"^  ®R"^ 
“R  or  “F 
non-dimensional 
non-dimensional 
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Machine  Program 

Program  Restrictions 

(1)  Thermal  capacitor  Indices  must  not  exceed  1000. 

(2)  There  must  be  no  more  than  20  CALL  AER0  statements. 

(3)  There  must  be  no  more  than  500  CALL  C0MC0N  statements. 

(4)  There  must  be  no  more  than  50  CALL  PL0T  statements. 

(5)  There  must  be  no  more  than  16  capacitors  In  any  CALL  PL0T  statement. 

Machine  Timing 

The  machine  time  may  generally  be  estimated  by  multiplying  the 
machine  time  per  time  step  by  the  number  of  time  steps  required  to  compute 
from  START  to  ST0P.  (See  the  last  part  of  Appendix  B  on  how  to  estimate 
the  number  of  time  steps.)  The  machine  time  per  time  step  may  be  estimated 
by  multiplying  the  time  required  by  each  Inner  loop  subroutine  by  the  number 
of  times  It  Is  called  and  summing  over  all  Inner  loop  subroutines. 

Subroutine  Name  Time 


F0RCER 

1 

mllll-second 

AER0 

20 

mllll-seconds 

RAD 

1 

mllll-second 

C^C0N 

5 

mllll-seconds 

CON 

1 

mllll-second 

CAP 

1 

mllll-second 

Machine  storage  Space  Requirements 

(a)  It  will  generally  prove  most  convenient  to  use  BCC  Library  Subroutines 
11.02.01-11.02.12  as  a  package.  In  the  somewhat  unlikely  event  that  space 
becomes  a  problem,  the  following  Information  may  be  useful: 


Number 

Name 

Locations 

11.02.01 

SET 

93 

11.02.02 

TRAJ 

40 

11.02.03 

AMBATM 

150 

11.02.04 

F0RCER 

40 

11.02.05 

AER0 

670 

11.02.06 

RAD 

85 

11.02.07 

C^C0N 

100 

11.02.08 

C0N 

65 

11.02.09 

CAP 

40 

11.02.10 

WRITE 

300 

11.02.11 

PL0T 

170 

11.02.12 

STEP 

140 

The  11.02.01-11.02.12  package  also  Includes  BCC  Library  Subroutines  DECIDE, 
PIFl  and  PIF2.  The  package  uses  6174  C0MC0N  locations. 
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(b)  Each  CALL  statement  In  the  control  program  requires  a  certain  number  of 
locations  depending  on  the  type  of  CALL  statement. 


Type  Locations 


CALL  SET 

5  + 

21 

CAI.L  TRAJ 

5 

CALL  AMBATM 

5 

CALL  F0RCER 

4 

CALL  AER0 

8  + 

31 

CALL  RAD 

6 

CALL  C0MC0N 

12 

CALL  C0N 

6 

CALL  CAP 

5 

CALL  WRITE 

2  + 

I 

CALL  PL0T 

4  + 

I 

CALL  STEP 

1 

IW  JOHM  HOKIM  UNWmity 

Amm  PMYiia  uioutmy 

Mva  nNo  lunuND 
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ADDENDUM 


This  addendum  Is  Included  In  order  to  instruct  the  reader  in  some 
parts  of  the  Applied  Physics  Laboratory  heat  transfer  program  which  need 
particular  care.  The  critical  aspect  of  the  Information  listed  below  has 
been  made  evident  to  the  BBE  heat  transfer  section  by  Its  experience  with 
numerous  computer  analyses  over  the  past  year.  The  problems  that  may  arise 
In  a  computer  program  run,  neglecting  obvious  technical  errors,  may  be 
broken  Into  three  categories: 

1.  Insufficient  Information  to  obtain  adequate  representation  of  data 
by  Interpolation.  This  problem  usually  arises  as  a  result  of  the 
fact  that  with  the  exception  of  the  AMBATM  temperature  function,  all 
tabular  functions  are  Interpolated  by  quadratic  equations. 

Example : 

Suppose  the  following  altitude  z  is  presented  as  a  function  of 
time  t. 


F2(l) 

=  0  . 

F2V(1)  -  N 

F2(2) 

=  2  . 

F2V(2)  =•  0 

F2(3) 

=•  10  . 

F2V(3)  =  1 

Zjj  F2(N)  =  100000. 

F2V(Nfl)  =*  100.  tj^ 

where  the  quadratic  equation 

2 

z  =  at  +  bt  +  c 

Is  used. 

At  t  -  tj  ■  0,  Zj^  ■  0  therefore  c  »  0 

at  t  ®  t^  “  1,  ^2  *  ^  therefore, 

2  =  a(l)^  +  b(l)  therefore, 

b  -  2  -  a 

at  t  -  t3  -  2,  10  -  a(2)2  +  (2 -a)  (2) 

5  -  2a  +  2-a 

therefore,  a  ■  3  and  b  ■  -  1, 

2 

and  z  ■  3t  -  t. 
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Now  to  find  the  minimum  value  of  z  In  the  Interval  t^  to  t^,  differentiate 
z  with  respect  to  t  and  set  equal  to  zero. 


dz 

dt 


therefore, 


and 


z  ■ 


0; 


1 

6 


since,  no  ambient  temperatures  and  pressures  are  listed  for 
negative  altitudes,  the  computer  will  give  faulty  results.  Had  more 
altitude  data  been  given  In  the  time  interval  from  t^  to  t^  then  no 
negative  z  values  would  occur  In  this  time  Interval.  Therefore,  the 
machine  run  would  be  good.  In  conclusion,  be  generous  with  data  at 
the  beginning,  at  Inflection  points  and  at  ends  of  functional  values. 
The  best  way  of  checking  out  a  heat  transfer  program  Is  to  allow  the 
computer  to  Initially  make  only  a  few  seconds  of  flight  time  analysis. 
That  Is  In  the  START  and  STOP  of  the  CALL  SET  statement ,  make  START 
0  seconds  and  make  STOP,  say,  2  seconds.  That  way,  If  something  Is 
wrong  In  the  Input  Information  corrections  can  be  made  without  wasting 
much  computer  time. 

2.  Inherent  limitations  of  convective  heat  transfer  equations.  The 

T'  -  Colburn  method  of  evaluating  the  coefficient  of  convective  heat 
transfer  Is  best  suited  to  free  stream  Mach  numbers  of  5  or  less. 

This  Is  the  method  built  into  the  Applied  Physics  Laboratory  program. 
Above  Mach  5  the  Van  Driest  method  of  evaluating  h  is  more  accurate. 
Also,  the  equations  for  h  must  be  modified  for  heat  transfer  at  stagna¬ 
tion  points.  This  Is  because  of  the  appearance  of  the  characteristic 
length  ^  In  the  denominator  of  the  equation  for  h. 

3.  Unwise  assumptions.  In  studies  Involving  plastic  surfaces.  It  must 
be  remeiid>ered  that  melting  and/or  ablating  may  take  place  at  high 
speeds.  The  APL  program  does  not  handle  the  changes  In  state  and 
geometry  that  may  thus  ensue.  Therefore,  computed  data  which  may 
look  good  on  surface,  may  be  useless  If  above  phenomena  had  taken 
place. 
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Another  condition  that  one  must  exercise  care  with  is  choosing 
adequately  thin  elements  for  the  representation  of  aerodynamic  surfaces 
made  of  fiberglass  phenolics  or  other  low  conductive  materials.  The 
temperature  Obtained  for  the  surface  element  is  used  as  the  wall  tempera¬ 
ture  in  the  computation  of  T*.  In  turn,  this  is  used  in  the  computation 
of  h.  For  the  case  of  low  conductive  materials,  when  the  surface  element 
is  thick,  then  will  be  smaller  during  the  ascent  of  a  missile  than  if 
the  surface  element  were  made  thin  because  of  the  larger  volume.  In  turn, 
the  coefficient  of  convective  heat  transfer  will  be  larger.  Thus,  over 
a  finite  missile  ascending  period  of  time,  more  heat  will  be  supplied  to 
the  thicker  surface  element.  Furthermore,  since  the  overall  thickness  of 
the  structure  which  is  receiving  aerodynamic  heating  is  a  fixed  dimension, 
the  inner  element  layers  will  eventually  have  to  absorb  the  excessive  heat 
supplied  to  the  surface  element.  Thus,  not  only  will  the  temperature  of 
the  surface  element  be  in  error,  but  also  the  temperatures  of  all  inner 
elements.  If  the  conductivity  of  the  material  is  large  as  in  the  case 
of  metals,  then  the  above  will  not  nearly  be  as  serious.  In  conclusion, 
for  temperature  gradient  supporting  materials,  generous  usage  of  break-up 
into  layers  should  be  made  in  the  mathematical  model. 
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